↜ Back to index Introduction to Numerical Analysis 1

Solution Lecture a4

Example solutions of Report 1.

Exercise 1

! Implementation of the midpoint method for y' = y sin t
implicit none
integer i, n
real y, h, t, f, k

n = 20
! interval (0, 5)
h = 5. / n

! Initial data
y = 1.
print *, 0., y

do i = 1, n
    ! One step of the midpoint method.
    t = (i - 1) * h
    f = y * sin(t)
    k = y + 0.5 * h * f
    t = (i - 0.5) * h
    f = k * sin(t)
    y = y + h * f
    print *, i * h, y
enddo
end

Exercise 2

! Implementation of the backward Euler method for y' = -10(y - cos t)
implicit none
integer i, n
real y, h, t, f

n = 10
h = 5. / n

! Initial data
y = 0.
print *, 0., y

do i = 1, n
    t = i * h
    y = (y + 10 * h * cos(t)) / (1 + 10 * h)
    print *, t, y
enddo

end

Exercise 3

1-2 without arrays

! Implementation of the midpoint method for a system y_1' = y_2, y_2' = - y_1
implicit none
integer i, n
real h, t
real y1, y2, k1, k2

n = 20
h = 10. / n

! Initial data
y1 = 1.
y2 = 0.
print *, 0., y1, y2

do i = 1, n
    k1 = y1 + 0.5 * h * y2
    k2 = y2 + 0.5 * h * (-y1)
    y1 = y1 + h * k2
    y2 = y2 + h * (-k1)

    print *, i * h, y1, y2
enddo

end

1-2 with arrays

! Implementation of the midpoint method for a system y_1' = y_2, y_2' = - y_1
implicit none
integer i, n
real h, t
real y(2), k(2)

n = 20
h = 10. / n

! Initial data
y = [1., 0.]
print *, 0., y

do i = 1, n
    k = y + 0.5 * h * [y(2), -y(1)]
    y = y + h * [k(2), -k(1)]

    print *, i * h, y
enddo

end

1-3

! Implementation of the midpoint method for the system y_1' = y_2, y_2' = - y_1
implicit none
integer i, n, m
real h, t
real y(2), k(2)

do m = 3, 10
    n = 2**m
    h = 10. / n

    ! Initial data
    y = [1., 0.]

    do i = 1, n
        k = y + 0.5 * h * [y(2), -y(1)]
        y = y + h * [k(2), -k(1)]

    enddo

    print *, n, sqrt((y(1) - cos(10.))**2 + (y(2) + sin(10.))**2)
enddo

end

Exercise 4

! Implementation of the backward Euler method for the system y_1' = y_2, y_2' = - y_1
implicit none
integer i, n
real h, t
real y(2)

n = 100
h = 10. / n

! Initial data
y = [1., 0.]
print *, 0., y

do i = 1, n
    y = (y + h * [y(2), -y(1)]) / (1 + h**2)

    print *, i * h, y
enddo

end